{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from urllib.request import urlretrieve\n",
    "import os\n",
    "\n",
    "import itk\n",
    "import numpy as np\n",
    "\n",
    "from itkwidgets import view, cm\n",
    "from ipywidgets import FloatProgress, Label, HBox, VBox, FloatText, ColorPicker, Button"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download input image\n",
    "file_name = 'Data/901-L.nrrd'\n",
    "if not os.path.exists(file_name):\n",
    "    url = 'https://data.kitware.com/api/v1/file/5ef373379014a6d84edb66f8/download'\n",
    "    urlretrieve(url, file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = itk.imread(file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56b8f1c6052d491fba4f8e49c30c7f6b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "camera = np.array([[-8.23076   , 16.761303  , 41.647186  ],\n",
    "       [ 4.855     ,  9.24      , 48.065998  ],\n",
    "       [ 0.3404956 , -0.19294305, -0.92023677]], dtype=np.float32)\n",
    "opacity_gaussians = [[{'position': 0.23888888888888893,\n",
    "   'height': 0.3090909090909091,\n",
    "   'width': 0.1722222222222221,\n",
    "   'xBias': 0.04708953460902432,\n",
    "   'yBias': 0}]]\n",
    "viewer = view(image,\n",
    "     cmap='bone_Matlab',\n",
    "     gradient_opacity=1.0,\n",
    "     background=(0.6,0.6,0.6),\n",
    "     ui_collapsed=True,\n",
    "     camera=camera,\n",
    "     opacity_gaussians=opacity_gaussians\n",
    "    )\n",
    "viewer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b0192ae2d00848fa934835b618650b37",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Segment bones'), FloatProgress(value=0.0, bar_style='info', description='progress:…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "label = Label('Segment bones')\n",
    "progress = FloatProgress(value=0.0,\n",
    "                         min=0.0,\n",
    "                         max=1.0,\n",
    "                         step=0.001,\n",
    "                         description='progress:',\n",
    "                         bar_style='info')\n",
    "progress\n",
    "box = HBox([label, progress])\n",
    "box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "itkPointSetD3 not loaded from module ITKRegistrationCommon because of exception:\n",
      " module 'itk.ITKRegistrationCommonPython' has no attribute 'itkPointSetD3'\n",
      "Template itk::LandmarkAtlasSegmentationFilter<itk::Image<signedshort,3>,itk::Image<signedshort,3>>\n",
      " already defined as <class 'itk.itkLandmarkAtlasSegmentationFilterPython.itkLandmarkAtlasSegmentationFilterISS3ISS3'>\n",
      " is redefined as <class 'itk.itkLandmarkAtlasSegmentationFilterPython.itkLandmarkAtlasSegmentationFilterISS3ISS3'>\n",
      "Warning: template already defined 'itk::LandmarkAtlasSegmentationFilter<itk::Image<signedshort,3>,itk::Image<signedshort,3>>'\n"
     ]
    }
   ],
   "source": [
    "bone_segmenter = itk.SegmentBonesInMicroCTFilter.New(image)\n",
    "cortical_thickness=0.1\n",
    "bone_segmenter.SetCorticalBoneThickness(cortical_thickness)\n",
    "\n",
    "def update_progress():\n",
    "    progress.value = bone_segmenter.GetProgress()\n",
    "bone_segmenter.AddObserver(itk.ProgressEvent(), update_progress)\n",
    "\n",
    "bone_segmenter.Update()\n",
    "\n",
    "progress.bar_style = 'success'\n",
    "progress.description = 'done.'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_image = bone_segmenter.GetOutput()\n",
    "label_file_name = file_name + '-labels.nrrd'\n",
    "itk.imwrite(label_image, label_file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9324dbdc9beb424b87b9d6ca7ab9ce85",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "label_image = bone_segmenter.GetOutput()\n",
    "view(image=image,\n",
    "     label_image=label_image,\n",
    "     opacity_gaussians=opacity_gaussians,\n",
    "     gradient_opacity=1.0,\n",
    "     background=(0.6,0.6,0.6),\n",
    "     camera=camera,\n",
    "     ui_collapsed=True,\n",
    "     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download atlas image\n",
    "atlas_file_name = '907-L.nrrd'\n",
    "if not os.path.exists(atlas_file_name):\n",
    "    url = 'https://data.kitware.com/api/v1/file/5ef375199014a6d84edb6a1f/download'\n",
    "    urlretrieve(url, atlas_file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlas = itk.imread(atlas_file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download atlas label image\n",
    "atlas_label_file_name = '907-L-atlas.nrrd'\n",
    "if not os.path.exists(atlas_label_file_name):\n",
    "    url = 'https://data.kitware.com/api/v1/file/5ef372559014a6d84edb627e/download'\n",
    "    urlretrieve(url, atlas_label_file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlas_label = itk.imread(atlas_label_file_name)\n",
    "# view(label_image=atlas_label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d5067c4b1054dc0b6526a5a599ea148",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "center_of_femur_head = np.array([4.072,11.347,43.816]).reshape((1,3))\n",
    "femur_shaft = np.array([7.713,9.107,42.649]).reshape((1,3))\n",
    "dent = np.array([3.787,11.172,45.173]).reshape((1,3))\n",
    "\n",
    "opacity_gaussians = [[{'position': 0.33611111111111114,\n",
    "   'height': 0.24545454545454548,\n",
    "   'width': 0.22499999999999998,\n",
    "   'xBias': 0.51,\n",
    "   'yBias': 0.4}]]\n",
    "view(atlas,\n",
    "     #label_image=atlas_label,\n",
    "     gradient_opacity=1.0,\n",
    "     background=(0.6,0.6,0.6),\n",
    "     cmap='bone_Matlab',\n",
    "     camera=camera,\n",
    "     label_image_blend=0.4,\n",
    "     opacity_gaussians=opacity_gaussians,\n",
    "     point_set_sizes=[15,]*3,\n",
    "     ui_collapsed=True,\n",
    "     point_set_representations=['spheres',]*3,\n",
    "     point_sets=[center_of_femur_head, femur_shaft, dent])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3bbac412e0b84989ba09c0653f32b7b7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "center_of_femur_head = np.array([4.945,9.225,46.011]).reshape((1,3))\n",
    "femur_shaft = np.array([7.424,9.474,44.711]).reshape((1,3))\n",
    "dent = np.array([5.018,9.058,47.535]).reshape((1,3))\n",
    "\n",
    "opacity_gaussians = [[{'position': 0.33611111111111114,\n",
    "   'height': 0.24545454545454548,\n",
    "   'width': 0.22499999999999998,\n",
    "   'xBias': 0.51,\n",
    "   'yBias': 0.4}]]\n",
    "viewer = view(image,\n",
    "     gradient_opacity=1.0,\n",
    "     background=(0.6,0.6,0.6),\n",
    "     cmap='bone_Matlab',\n",
    "     camera=camera,\n",
    "     opacity_gaussians=opacity_gaussians,\n",
    "     point_set_sizes=[15,]*3,\n",
    "     point_set_representations=['spheres',]*3,\n",
    "     point_sets=[center_of_femur_head, femur_shaft, dent])\n",
    "viewer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0fb7e154557b4ed68bac29748b3b0a19",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(VBox(children=(HBox(children=(ColorPicker(value='#d60000', concise=True, description='Head Cent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head_center_widgets_top = [ColorPicker(description='Head Center ', concise=True, value='#d60000', disabled=True),\n",
    "                           Button(description='From last click', tooltip='Set from last clicked slice point', icon='dot-circle-o'),]\n",
    "head_center_widgets_bottom = [FloatText(value=4.945, description='x:'),\n",
    "                              FloatText(value=9.225, description='y:'),\n",
    "                              FloatText(value=46.011, description='z:'),]\n",
    "def head_center_from_click(change):\n",
    "    position = viewer.clicked_slice_point.position\n",
    "    for i in range(3):\n",
    "        head_center_widgets_bottom[i].value = position[i]\n",
    "head_center_widgets_top[1].on_click(head_center_from_click)\n",
    "head_center_box = VBox([HBox(head_center_widgets_top), HBox(head_center_widgets_bottom)])\n",
    "\n",
    "shaft_widgets_top = [ColorPicker(description='Shaft', concise=True, value='#8c39ff', disabled=True),\n",
    "                     Button(description='From last click', tooltip='Set from last clicked slice point', icon='dot-circle-o'),]\n",
    "shaft_widgets_bottom = [FloatText(value=7.424, description='x:'),\n",
    "                        FloatText(value=9.474, description='y:'),\n",
    "                        FloatText(value=44.711, description='z:'),]\n",
    "shaft_box = VBox([HBox(shaft_widgets_top), HBox(shaft_widgets_bottom)])\n",
    "\n",
    "indent_widgets_top = [ColorPicker(description='Indent', concise=True, value='#018700', disabled=True),\n",
    "                     Button(description='From last click', tooltip='Set from last clicked slice point', icon='dot-circle-o'),]\n",
    "indent_widgets_bottom = [FloatText(value=5.018, description='x:'),\n",
    "                         FloatText(value=9.058, description='y:'),\n",
    "                         FloatText(value=47.535, description='z:'),]\n",
    "indent_box = VBox([HBox(indent_widgets_top), HBox(indent_widgets_bottom)])\n",
    "\n",
    "def update_positions():\n",
    "    head_center_point = np.array([ft.value for ft in head_center_widgets_bottom]).reshape((1,3))\n",
    "    shaft_point = np.array([ft.value for ft in shaft_widgets_bottom]).reshape((1,3))\n",
    "    indent_point = np.array([ft.value for ft in indent_widgets_bottom]).reshape((1,3))\n",
    "    viewer.point_sets = [head_center_point, shaft_point, indent_point]\n",
    "\n",
    "\n",
    "def head_center_from_click(change):\n",
    "    if viewer.clicked_slice_point is None:\n",
    "        return\n",
    "    position = viewer.clicked_slice_point.position\n",
    "    for i in range(3):\n",
    "        head_center_widgets_bottom[i].value = position[i]\n",
    "    update_positions()\n",
    "head_center_widgets_top[1].on_click(head_center_from_click)\n",
    "def shaft_from_click(change):\n",
    "    if viewer.clicked_slice_point is None:\n",
    "        return\n",
    "    position = viewer.clicked_slice_point.position\n",
    "    for i in range(3):\n",
    "        shaft_widgets_bottom[i].value = position[i]\n",
    "    update_positions()\n",
    "shaft_widgets_top[1].on_click(shaft_from_click)\n",
    "def indent_from_click(change):\n",
    "    if viewer.clicked_slice_point is None:\n",
    "        return\n",
    "    position = viewer.clicked_slice_point.position\n",
    "    for i in range(3):\n",
    "        indent_widgets_bottom[i].value = position[i]\n",
    "    update_positions()\n",
    "shaft_widgets_top[1].on_click(shaft_from_click)\n",
    "for float_input in head_center_widgets_bottom + shaft_widgets_bottom + indent_widgets_bottom:\n",
    "    float_input.observe(lambda x: update_positions(), 'value')\n",
    "\n",
    "position_widgets = VBox([head_center_box, shaft_box, indent_box])\n",
    "position_widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "landmark_registration = itk.LandmarkAtlasSegmentationFilter[type(image), type(label_image)].New()\n",
    "landmark_registration.SetInput(image)\n",
    "landmark_registration.SetInput(1, atlas)\n",
    "landmark_registration.SetInputLabels(label_image)\n",
    "landmark_registration.SetAtlasLabels(atlas_label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "LandmarksType = itk.vector[itk.Point[itk.D, 3]]\n",
    "\n",
    "input_landmarks = LandmarksType()\n",
    "input_landmarks.push_back([ft.value for ft in head_center_widgets_bottom])\n",
    "input_landmarks.push_back([ft.value for ft in shaft_widgets_bottom])\n",
    "input_landmarks.push_back([ft.value for ft in indent_widgets_bottom])\n",
    "landmark_registration.SetInputLandmarks(input_landmarks)\n",
    "\n",
    "atlas_landmarks = LandmarksType()\n",
    "atlas_landmarks.push_back([4.072,11.347,43.816])\n",
    "atlas_landmarks.push_back([7.713,9.107,42.649])\n",
    "atlas_landmarks.push_back([3.787,11.172,45.173])\n",
    "landmark_registration.SetAtlasLandmarks(atlas_landmarks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 44.6 s, sys: 6.34 s, total: 50.9 s\n",
      "Wall time: 12.6 s\n"
     ]
    }
   ],
   "source": [
    "%time landmark_registration.Update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cb660df12ce849f0af5e2b4f86431674",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(image=image,\n",
    "     cmap=cm.bone,\n",
    "     camera=camera,\n",
    "     gradient_opacity=1.0,\n",
    "     background=(0.6,0.6,0.6),\n",
    "     opacity_gaussians=opacity_gaussians,\n",
    "     label_image=landmark_registration.GetOutput())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_masked_trabecular = itk.mask_negated_image_filter(image, mask_image=label_image, masking_value=2)\n",
    "input_masked_cancellous = itk.mask_negated_image_filter(image, mask_image=label_image, masking_value=1)\n",
    "input_masked = itk.add_image_filter(input_masked_trabecular, input_masked_cancellous)\n",
    "del input_masked_trabecular\n",
    "del input_masked_cancellous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlas_masked = itk.mask_image_filter(atlas, mask_image=atlas_label, masking_value=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_masked_trabecular = itk.mask_negated_image_filter(label_image, mask_image=label_image, masking_value=2)\n",
    "label_masked_cancellous = itk.mask_negated_image_filter(label_image, mask_image=label_image, masking_value=1)\n",
    "label_masked = itk.add_image_filter(label_masked_trabecular, label_masked_cancellous)\n",
    "label_masked = itk.cast_image_filter(label_masked, ttype=(type(label_masked), itk.Image[itk.UC, 3]))\n",
    "del label_masked_trabecular\n",
    "del label_masked_cancellous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_masked_map = itk.label_image_to_label_map_filter(label_masked)\n",
    "cropped = itk.auto_crop_label_map_filter(label_masked_map, crop_border=[6,]*3)\n",
    "label_masked_cropped = itk.label_map_to_label_image_filter(cropped)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e5979f5403a6486da2d1eed120742620",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "input_masked_cropped = itk.extract_image_filter(input_masked,\n",
    "                                                extraction_region=label_masked_cropped.GetLargestPossibleRegion())\n",
    "view(input_masked_cropped,\n",
    "     label_image=label_masked_cropped,\n",
    "     cmap=cm.bone,\n",
    "     camera=camera,\n",
    "     gradient_opacity=1.0,\n",
    "     background=(0.6,0.6,0.6),\n",
    "     opacity_gaussians=opacity_gaussians\n",
    "     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = landmark_registration.GetFinalTransform()\n",
    "atlas_affine_transformed = itk.resample_image_filter(atlas_masked,\n",
    "                                                     use_reference_image=True,\n",
    "                                                     reference_image=input_masked_cropped,\n",
    "                                                     transform=transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "nearest_interpolator = itk.NearestNeighborInterpolateImageFunction.New(atlas_label)\n",
    "atlas_labels_affine_transformed = itk.resample_image_filter(atlas_label,\n",
    "                                                            use_reference_image=True,\n",
    "                                                            reference_image=input_masked_cropped,\n",
    "                                                            transform=transform,\n",
    "                                                            interpolator=nearest_interpolator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "563407adb73540d08b178f418f19a2d4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "itk.imwrite(input_masked_cropped, file_name + '-input-masked-cropped.nrrd')\n",
    "itk.imwrite(atlas_affine_transformed, file_name + '-atlas-affine-transformed.nrrd')\n",
    "itk.imwrite(atlas_labels_affine_transformed, file_name + '-atlas-labels-affine-transformed.nrrd')\n",
    "view(atlas_affine_transformed,\n",
    "     label_image=atlas_labels_affine_transformed,\n",
    "     cmap=cm.bone,\n",
    "     camera=camera,\n",
    "     gradient_opacity=1.0,\n",
    "     background=(0.6,0.6,0.6),\n",
    "     opacity_gaussians=opacity_gaussians\n",
    "     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_masked_cropped_f = itk.imread(file_name + '-input-masked-cropped.nrrd', itk.F)\n",
    "atlas_affine_transformed_f = itk.imread(file_name + '-atlas-affine-transformed.nrrd', itk.F)\n",
    "atlas_labels_affine_transformed_f = itk.imread(file_name + '-atlas-labels-affine-transformed.nrrd', itk.F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3min 37s, sys: 5.79 s, total: 3min 43s\n",
      "Wall time: 49.3 s\n"
     ]
    }
   ],
   "source": [
    "%time registered_atlas, transform_parameter_object = itk.elastix_registration_method(input_masked_cropped_f, atlas_affine_transformed_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "transformer = itk.TransformixFilter[type(atlas_labels_affine_transformed_f)].New()\n",
    "transformer.SetMovingImage(atlas_labels_affine_transformed_f)\n",
    "transformer.SetTransformParameterObject(transform_parameter_object)\n",
    "transformer.Update()\n",
    "transformed = transformer.GetOutput()\n",
    "transformed_cast = itk.cast_image_filter(transformed, ttype=(type(transformed), itk.Image[itk.UC, 3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "34e7066ba7b043c2842c4f077e49772e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(background=(0.6, 0.6, 0.6), camera=array([[-8.23076   , 16.761303  , 41.647186  ],\n",
       "       [ 4.855     ,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "opacity_gaussians = [[{'position': 0.5,\n",
    "   'height': 0.5636363636363637,\n",
    "   'width': 0.5,\n",
    "   'xBias': 0.026666666666666672,\n",
    "   'yBias': 0}]]\n",
    "registered = view(input_masked_cropped_f,\n",
    "                  label_image=transformed_cast,\n",
    "                  label_image_blend=0.7,\n",
    "                  camera=camera,\n",
    "                  gradient_opacity=1.0,\n",
    "                  background=(0.6,0.6,0.6),\n",
    "                  cmap=cm.bone,\n",
    "                  opacity_gaussians=opacity_gaussians,\n",
    "                  ui_collapsed=True,\n",
    "                  )\n",
    "registered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = np.ones((10,), dtype=np.float32)*0.1\n",
    "weights[5] = 1.0\n",
    "registered.label_image_weights = weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "roi = itk.mask_negated_image_filter(transformed_cast, mask_image=transformed_cast, masking_value=5)\n",
    "# view(label_image=roi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "morphometry_filter = itk.BoneMorphometryFeaturesFilter.New(input_masked_cropped_f)\n",
    "morphometry_filter.SetMaskImage(roi)\n",
    "morphometry_filter.Update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BVTV 0.7600745916651215\n",
      "TbN 15.203703884426831\n",
      "TbTh 0.049992725288715116\n",
      "TbSp 0.015780720945284547\n",
      "BSBV 40.00582061589391\n"
     ]
    }
   ],
   "source": [
    "print('BVTV', morphometry_filter.GetBVTV())\n",
    "print('TbN', morphometry_filter.GetTbN())\n",
    "print('TbTh', morphometry_filter.GetTbTh())\n",
    "print('TbSp', morphometry_filter.GetTbSp())\n",
    "print('BSBV', morphometry_filter.GetBSBV())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
